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Abstract 

Past molecular dynamics studies of thermal transport have predominantly used Stillinger- Weber 
potentials. As materials continuously shrink, their properties increasingly depend on defect and 
surface effects. Unfortunately, Stillinger- Weber potentials are best used for diamond-cubic-like bulk 
crystals. They cannot represent the energies of many metastable phases, nor can they accurately 
predict the energetics of defective and surface regions. To study nanostructured materials, where 
these regions can dominate thermal transport, the accuracy of Tersoff potentials in representing 
these structures is more desirable. Based upon an analysis of thermal transport in a GaN system, 
we demonstrate that the cutoff function of the existing Tersoff potentials may lead to problems in 
determining the thermal conductivity. To remedy this issue, improved cutoff schemes are proposed 
and evaluated. 
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Thermal conductivity of semiconductors is often the limiting factor for important techno- 
logical applications. For instance, a higher thermal conductivity improves heat dissipation, 
which also enables an increase in the microelectronic device density [T] . On the other-hand, 
a lower thermal conductivity has long been sought in order to push the energy conversion ef- 
ficiency of thermoelectric devices beyond a critical high value |2J. As the system dimensions 
continuously shrink towards the nanoscale, the thermal conductivity increasingly depend on 
defects and surfaces of materials. Because experimental study is challenging at the nanoscale, 
a theoretical understanding of the influence of defects and surfaces on thermal conductivity 
is essential for optimizing nanoscale applications. With the structures and energetics of de- 
fects and surfaces precisely captured in a computational lattice, molecular dynamics (MD) 
has emerged as a powerful means to study the thermal transport in nano crystals. MD has 
been successfully applied for both covalent (semiconductors) [3HT6] and other (e.g., ionic) 
[T71EH] systems. 

MD simulations of thermal transport in semiconductors have used predominantly the 
Stillinger- Weber [25] rather than the Tersoff/ Brenner [2SH28] types of potentials. Of the 
papers impartially cited here, a majority [3HT3] applied Stillinger- Weber potentials and only 
a few [HHT6] used Tersoff/Brenner types of potentials. Recently, we published a series 
of papers on MD simulations of thermal transport in GaN [T0T[T2] . While these studies 
employed a Stillinger- Weber GaN potential [29, 30J, significant efforts were also made to 
explore the Tersoff GaN potential [31] during the course of the work. With this particular 
parameterization of the Tersoff potential, we found it difficult to ensure energy conservation 
during a "constant energy" MD simulation unless the time step size was significantly reduced 
from the 0.001 ps we normally used for other potentials, and the thermal conductivities 
obtained were more than one order of magnitude smaller than the values we would otherwise 
obtain from the Stillinger- Weber potential [291 E0] • 

From a physics point of view, Stillinger- Weber potentials are empirically designed for 
diamond-cubic-like bulk environment and they do not work well in defective and surface 
regions. For instance, the angular term in Stillinger- Weber potential is essentially a positive 
parabolic energy penalty for deviations of bond angles from the tetrahedral angle. As this 
positive energy term is repulsive, any structures with non-tetrahedral bond angles would 
have longer bond lengths than those in diamond-cubic-like crystals or dimers. This is not 
always true in reality. Furthermore, under the condition that the diamond-cubic-like crystal 
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has the lowest energy, the angular term in Stillinger- Weber potential must be constrained 
above a certain level. As a result the potential cannot correctly describe the trend of 
energies in many other metastable phases. Tersoff potentials, on the other hand, inherit the 
quantum mechanics-based bond order concept [32ti3l] . They are therefore well suited for 
treating defects and surfaces. In particular, we note that the GaN Tersoff potential [31] was 
carefully parameterized using cohesive energies, atomic volumes, defect properties, melting 
temperatures, and solubility limits of a variety of stable and metastable phases. Its fidelity 
has also been validated in MD simulations of vapor deposition that correctly showed the 
crystalline growth of the lowest energy phase [35J . 

The energy non-conservation and the abnormally low thermal conductivity are not trans- 
parent problems. Intuitively, they cannot be attributed directly to the potential format. 
Resolving these significant problems is important to the materials modeling community, as 
it will both promote the use of Tersoff potential in future MD simulations and improve the 
atomistic simulations of semiconductor systems. Understanding the abnormally low thermal 
conductivity predicted by the GaN Tersoff [3T] potential may also result in the discovery of 
a new thermal scattering phenomenon that can be utilized to engineer thermal conductivity 
of materials. This paper will address these issues. 

Following the format of the GaN paper [31], the Tersoff potential is expressed as 

E = J2Y1 ^ ~ Bij ■ fai ( r v)] (!) 

i j>i 

where the bond order is defined as 

with the local angular-dependent variable Xij calculated as 

Xij = X] $ ik ( r * fc ) ' 9ik (^' ifc ) ' exp t 2/iiA: ( r ^' ~~ Tik ^ fi) 

Eqs. Q and ^ involve three pair functions 0^, ipij, and representing respectively the 
i-j repulsion, the i-j attraction, and the i-k cutoff (in the angularly dependent x term). Here 
the symbol represents the usual relative distance between atoms i and j and 0^ is the 
angle between the triplet jik with i at its vertex. The cutoff function / is taken to be a 
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sinoidal spline: 
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^ij ^ ^ij — Rij ^ij 
0, T^j Rij ~\~ -^ij 

Here r s = R — D is the cutoff starting radius and r c = R + D is the cutoff distance at which 
point the interactions go to zero. To impose the cutoff of the <fi and ip functions, the Tersoff 
potentials simply define <f>ij(r r .j) = ./',;(/•;,) • V$ (r^-) and ipij{r rd ) = ■ V^ry), where 

V^r) and V (r) are two exponentially decaying functions |26ti28l 131] that are not subject 
to the discussions here. Due to the use of Eq. Q, the three functions 0, ip, and / have 
continuous values and the first derivatives, but their second derivatives abruptly change at 
the junction points r s and r c . This is contrary to the Stillinger- Weber potentials whose high 
order derivatives are always continuous. 

By dividing the r s - r c range evenly into three sections with two internal points 
ri,r 2 G (r s ,r c ), a cubic spline cutoff approach can be used to remove the discontinuous 
2nd-derivatives of the 0, ip, and / functions. Using <ft as an example, we can write, 



^ij ^ ^s,ij 



En=2 Vff ( r s,ij) / \n / \3 

n=0 n\ V *J ~~ rs >*i/ "i fl 3 Vij ~ r s,ij) j r s,ij < f"ij < T\^j 



V n=3 6 (r ■■ 

Z^n=0 u n V %3 
c 3 ( r ij ~ r c,ij, 

0. 



T %ij) , 



n,ij < rij < r 2 - 



T"ij ^* '™C,ij 



(5) 



Eq. (|5]) is continuous up to the second derivative at r s and r c . It involves six parameters a 3 , 
&o - &3 and C3 that can be solved by requiring the values as well as the first and the second 
derivatives to be continuous at the two added junction points 7*1 and r 2 . For an example 
comparison, the modified and unmodified 4>GaGa — ipGaGa and fcaGa functions are shown as 
the black and red dash lines respectively in Fig. [TJ Note that equating the bond energy 
with 4>GaGa — ipGaGa implies a bond order of B GaGa = 1. This does not represent the real 
scenario but is sufficient for the discussion here. It can be seen that the modified potential 
does not significantly differ from the original potential. In fact, Eq. ^ should be applicable 
to any existing Tersoff potentials without a need to re-parameterize. This cubic spline cutoff 
scheme has been implemented in the MD simulator LAMMPS [36] . 
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FIG. 1: GaGa pair interaction functions (4>GaGa - 4>GaGa and fcaGaj- 

We suspect that the discontinuous 2nd-derivatives are likely to be the source of the 
problems with the Tersoff potentials. To explore this, we calculate the bond energy <pij (r^) — 
Bij ■ ipij as a function of bond length for different bonds ij = GaN, GaGa, NN using a 
fixed bond order estimated from the equilibrium wurtzite GaN and Eq. Q . The results 
are shown in Fig. [2] using black, long red dash, and short blue dash curves to represent, 
respectively, the Ga-N, Ga-Ga, and N-N interactions. The radial distribution of atoms was 
also determined from a GaN wurtzite equilibrated at a temperature of 300 K, and the result 
is included in Fig. [2] as the shaded area. Fig. [2] indicates that the third- or higher- nearest 
neighbors do not contribute to material properties as the corresponding radial distributions 
do not fall within ranges of any of the three interaction curves. Similarly, the N-N interaction 
does not contribute to properties because the N-N energy curve only overlaps with the 
nearest neighbor radial distribution which is between Ga and N. In fact, only the nearest 
neighbor Ga-N and the second-nearest neighbor Ga-Ga interactions contribute to properties. 
However, because the cutoff starting point r s ^ a N (= 2.7 A [31]) is well above the upper limit 
of the first radial distribution peak, how the Ga-N interaction is truncated has no effect on 
properties. Clearly, the effect due to a discontinuous 2nd- derivative, if any, must come from 
the Ga-Ga interaction. 

In addition to the cubic spline cutoff modification, Fig. [2] suggests numerous other ways 
to remove or mitigate the discontinuous 2nd-derivatives of the functions. As indicated 
by the circle in the figure, the Ga-Ga interaction has an insignificant contribution to the 
energy near the cutoff distance. In fact, there is no energy contribution at all at K 
because the second neighbor (Ga-Ga) distribution would then reduce to a single line at 
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FIG. 2: Different pair interactions in wurtzite GaN crystal. 

approximately the peak position outside of the Ga-Ga interaction range. This means that 
neglecting the Ga-Ga interactions would only cause a negligible energy change and would 
not alter the relative phase stabilities. Similarly, a reduction of the cutoff starting point 
r s ,GaGa (at a constant r c ^aGa cutoff distance) also does not significantly change the energy. 
A reduced r s ^GaGa, however, increases the cutoff transition range and therefore can be used to 
reduce magnitude of the 2nd- derivative near the transition. When no Ga-Ga interactions are 
considered, the only remaining interactions are between Ga and N, whose radial distribution 
is near the Ga-N interaction energy minimum which is far away from the cutoff distance. 
Consequently, the Ga-N neighbors do not change during thermal transport simulations where 
no phase transformation occurs. Hence, we can initialize the Ga-N neighbors at the start 
of simulations and not perform reneighboring calculations during the rest of simulations. It 
should be noted that the purpose of cutoff in atomistic simulations is to define the neighbor 
lists. When the correct Ga-N neighbors are initialized and are not recalculated, there is 
really no need to apply the cutoff algorithm to the Ga-N interaction functions. Hence, we can 
simply set fcaN = 1 to remove the discontinuous 2nd-derivatives of the Ga-N interactions. 

Based upon these observations, several schemes (as illustrated in Table [j]) are designed 
to test the effects of the discontinuous 2nd-derivatives of the Ga-Ga functions near their 
cutoffs: A, the original potential without modifications (time step size 0.0010 ps unless 
indicated otherwise); B, the original potential simulated at a reduced time step size of 0.0005 
ps; C, a reduced r s GaGa value to reduce the 2nd-derivatives (and hence the extent of the 
discontinuity) of the Ga-Ga pair functions at r c ^GaGa'i D, with only Ga-N interactions that 
are initialized without reneighboring; E, same as D with a further condition of fcaN = 1, 
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TABLE I: Results of temperature drift AT (K) and thermal conductivity k (W/K ■ m) obtained 
from MD simulations over a 10 ns period using different combinations of GaGa cutoff function, 



GaGa cutoff starting point r. 




reneighboring 


and timp stpn 


S1ZPS dt I T)S ) 




schemes 


A 


13 


c 


D 


ili 


1 


cutoff function 


sinoidal 


sinoidal 


sinoidal 


sinoidal 


none (f = 1) 


cubic spline 


r s ,GaGa 


2.72 


2.72 


2.32 


2.72 


N/A 


2.72 


reneighboring 


yes 


yes 


yes 


no 


no 


yes 


dt 


0.0010 


0.0005 


0.0010 


0.0010 


0.0010 


0.0010 


AT 


27 


3 


23 


3 


<1 


<1 


K 


5.4±0.1 


6.2±0.04 


19.0±0.3 


13.9+0.1 


14.6±0.1 


12.1±0.1 



and F, the cubic spline modification of the potential with other conditions the same as in 
scheme A. All these schemes were explored in thermal transport simulations. 

The "direct method" MD approach [T0] - [T2] was employed to calculate the thermal con- 
ductivities of GaN wurtzite crystal along the [0001] direction. The computational system 
used for the simulations is shown in Fig. [3j where the red and blue colors indicate respec- 
tively the hot and cold regions. The systems are 136, 50, and 6 cells in the x-, y-, and z- 
directions, corresponding approximately to a 708x276x18 A 3 material volume. Four voids 
are evenly created between the hot and the cold regions by removing six Ga and six N atoms 
as shown in Fig. |3j All simulations were carried out at a temperature of 300 K using peri- 
odic boundary conditions in all three coordinate directions. The simulated system therefore 
represents an infinite GaN crystal containing uniform defects. Thermal conductivities were 
averaged over 10 ns. The resulting temperature drift AT (energy drift) and thermal con- 
ductivity are compared as a function of simulation conditions in Table [I] for various schemes. 

Table [J indicates that all simulations that involve discontinuous 2nd-derivatives of the 
Ga-Ga interactions, such as schemes A and C, have an energy conservation problem at 
a time step size of 0.0010 ps (the temperature drift during the 10 ns simulation is AT 
> 20 K). This finding agrees well with the previous result [37] that non-smoothness of the 
potential functions causes energy drift. Reducing the time step to 0.0005 ps can significantly 
improve the energy conservation, as in scheme B. The thermal conductivities obtained from 
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FIG. 3: Atomistic configuration for MD simulations. 

simulations with large energy drift may not be accurate. Nonetheless, Table [I] clearly shows 
that the finite 2nd-derivatives of the Ga-Ga functions at the cutoff distance also contribute 
to the abnormally low thermal conductivities, as with schemes A and B. Reducing the 2nd- 
derivatives by reducing the cutoff starting point r St G a Ga increases thermal conductivity, as 
in scheme C. All simulations that remove the discontinuous 2nd-derivatives, schemes D - E, 
have good energy conservation and reasonably large thermal conductivities. Note that both 
schemes D and E ignore the Ga-Ga interactions and use the Ga-N neighbors determined at 
the start of simulations. Scheme E further assumes fcaNijcaN) = 1- However, this change 
of the Ga-N cutoff does not have a significant effect because the Ga-N spacing of the crystal 
is well below the cutoff starting point r s ,GaN as mentioned above. As a result, schemes D and 
E are expected to produce similar results. Because the reduction in thermal conductivity 
due to non-smooth cutoff function is an unphysical artifact, and because scheme E does not 
apply any cutoff modifications, the good agreement between the results of schemes E and F 
strongly validates the cubic spline cutoff approach. 

In summary, we discovered that the existing Tersoff potential formats do not always 
produce the correct results due to the discontinuous 2nd-derivatives of the cutoff functions. 
This problem can be solved by applying the cubic spline cutoff modifications. In some 
thermal transport simulations where atomic neighbors do not change, one can initialize the 
neighbor lists at the start of simulations and not perform reneighboring. This is not only 
computationally efficient, but also eliminates the need to apply the cutoff algorithms to the 
functions and thereby removes the problem of the Tersoff potential. When the discontinuous 
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2nd-derivatives come from an insignificant second neighbor interaction, the problem can be 
simply solved by ignoring that interaction. Finally, we discovered the fundamental fact 
that the finite 2nd-derivatives of the interaction between atoms near the cutoff distance can 
significantly change the thermal conductivity of the lattice. This observation calls further 
studies of the fidelity of interatomic potentials where the cutoff distance is often chosen rather 
than determined from fundamental insights. Ab initio studies to remedy this issue should 
enable molecular dynamics to reliably predict properties, such as thermal conductivity, at 
the nanoscale and promote new material discovery. 

The discontinuous second-derivatives lead to jumps in elastic constants. This would pre- 
vent realistic studies of mechanical properties of semiconductors (in particular, the plastic 
deformation and the dislocation behavior). The proposed modification of the Tersoff poten- 
tials can therefore impact many problems including the mechanical properties of materials. 

Sandia is a multi-program laboratory operated by Sandia Corporation, a Lockheed Martin 
Company, for the United States Department of Energy National Nuclear Security Adminis- 
tration under contract DEAC04-94AL85000. This work was performed under a Laboratory 
Directed Research and Development (LDRD) project. We are also grateful for helpful dis- 
cussions with A. P. Thompson and G. J. Wagner. 
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